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Abstract 

This extended abstract further develops the algorithms in [8] and 
[9] for rewriting expressions involving differential operators. The dif- 
ferential operators that we have in mind arise in the local analysis of 
nonlinear dynamical systems. In this work, we extend these algorithms 
in two different directions: YVe generalize the algorithms so that they 
apply to differential operators on groups and we develop the data struc- 
tures and algorithms to compute symbolically the action of differential 
operators on functions. Both of these generalizations are needed for 
applications. 
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Abstract 

In this paper, we review some recent work relating families of trees 
to symbolic algorithms for the exact computation of series which ap- 
proximate solutions of ordinary differential equations. It turns out that 
the vector space whose basis is the set of finite, rooted trees carries a 
natural multiplication related to the composition of differential opera- 
tors, making the space of trees an algebra. This algebraic structure can 
be exploited to yield a variety of algorithms for manipulating vector 
fields and the series and algebras they generate. 

1 Introduction 

In this paper, we review some recent work relating families of trees to sym- 
bolic algorithms for the exact computation of series which approximate so- 
lutions of ordinary differential equations. It turns out that the vector space 
whose basis is the set of finite, rooted trees carries a natural multiplication 
related to the composition of differential operators, making the space of trees 
an algebra. This algebraic structure can be exploited to yield a variety of 
algorithms for manipulating vector fields and the series and algebras they 
generate. t 

In Section 3, we introduce and explore the algebraic structure of trees. 
Section 4 describes a simplification algorithm for the rewriting of symbolic 
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expressions involving vector fields. Section 5 describes an algorithm for 
generating explicitly tntegrable flows associated with nilpotent Lie algebras. 
Section 6 exploits the relation between Taylor series and trees to study a 
class of intrinsic numerical integrators. We begin in Section 2 with some 
background. 

The results surveyed here are the work of a variety of mathematicians: 

[ especially want to mention the contributions of my collaborators Peter 
Crouch, Matthew Grayson and Richard Larson. The work on algebras of 
trees and its applications to symbolic computation is joint work with Richard 
Larson. All of the algorithms described here rest upon this foundation. The 
work on explicitly integrable flows and nilpotent Lie algebras is joint work 
with Matthew Grayson. The work on numerical algorithms evolving on 
groups is joint work with Peter Crouch. 

2 Background 
Consider a differential equation 

*(0 = £i(*(0) + «£j(*(0). x( 0) = x°eR N , (i) 

where E \ and Ej are vector fields and u is a parameter. In applications, 
u will be either a small perturbation u = c, a control t — ► or simply 

the constant u = 1. Unless the vector fields E\ and £2 are very special, no 
algorithm is known which will return the general solution to the system in 
closed form. Our objective is to find efficient algorithms to compute various 
approximate solutions of the differential equation exactly using symbolic 
computation. 

Although the impact of symbolic computation in this area is recent, the 
connection between the existence of closed form solutions and the approxi- 
mation of general solutions is a traditional theme, dating back to at least the 
the nineteenth century. One can distinguish two approaches. One, cham- 
pioned by Lie, is based upon algebra and geometry and concerns us here; 
the other, championed by Weirstrass and Poincare, is based upon complex 
function theory. * 

Consider a group of transformations acting on of the form 

• j 1 > ♦ • • 1 ^r)» — T • • • » 

with the property that the group permutes the solutions of the nonlinear 
system (1). Lie asked the question [51] and [52], How can information about 
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the tranformation group be used to help integrate the differential equation i 
To answer this question, Lie introduced the infinitesimal generators of the 
group 


M*) 


ds k dx» 


1 < k < r 


and showed that the Ak satisfy 


[A % ,Aj] - ^2 c xjAk> 
i 


where [■,■] is the commutator, or Lie bracket , 


[Ai,A } ] = AiAj - A,A i% 

and the cf are constants. For example, Lie showed that if there is a one 
parameter group of transformations permuting the solutions of a nonlinear 
system in the plane (x,f), then the integrating factor for the equation may 
be read off from the infinitesimal generator. 

Since Lie's time, this basic question has contributed to the development 
of a number of different fields: 


• The vector fields ,4 ; generate a filtered Lie algebra, which is usually 
infinite dimensional, and is the infinitesimal version of the continu- 
ous pseudogroup of transformations generated by the <J>,. Prior work 
has focused on the geometry and structure theory of these algebras; 
important contributions have been made by Guillemin and Sternberg 
[37] and [36], and Hermann [52], [53], building upon the earlier work 
of Cartan, Ehresmann and Spencer. 

• Formal sums of iterated powers of vector fields, or Lie series, have 
been developed by Grobner [25], [26] and Knapp and Wanner [46], 
[47] into an operational calculus and used to approximate the solu- 
tions of differential equations. Lie transform methods have also been 
used in perturbation theory by Rand [59] and Meyer [54], in celestial 
mechanics by Deprit [18}, and in particle physics by Dragt [20]. 

• Explicit series computations of solutions of differential equations have 
a number of interesting connections with combinatorics. Chen [13] 
makes uses of the shuffle product, Joyal [45], Labelle [49], Leroux [50], 
and Viennot [71] employ trees and species, while Rota, Krahaner and 
Odlyzko [62] exploit the umbral calculus. 
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• Ritt and Kolchin, building upon earlier work of Picard and Vessiot, 
developed the field of differential Galois theory. The goal is to obtain a 
theory describing the solvability of differential equations analogous to 
Galois’ theory describing the solvability of algebraic equations. The 
survey by Singer [67] provides a good description of this field from 
the viewpoint of symbolic computation. Other relevant contributions 
include the beginnings of a differerential Grobner theory [3], [38], anal- 
ogous to the Grobner theory in commutative algebra; and a Hopf al- 
gebraic interpretation of Picard- Vessiot theory by Takeuchi [70]. 

• There is now a resurgence of interest in using symmetry groups to 
help integrate differential equations. This direction of research has 
been active in the Soviet Union for some time, but during the past 
decade there has been increased interest in the United States. Impor- 
tant contributions have been made by Olver [57], Schwarz [64], and 
Dluman and Cole [5]. Closely related is the study by Caviness [66] of 
conservation laws for differential equations. 

During the past several years, there has been increasing interest in sym- 
bolic computation and differential equations. Work has proceeded in a num- 
ber of directions, and is based upon both the Lie and the Weifstrass and 
Poincare traditions: 

• Zippel [73] is writing a modular symbolic computation system which 
supports the ability to call high quality numerical routines. Using such 
a system, he has shown how symbolic algorithms can be used to select 
appropriate numerical algorithms. 

• Guckenhemier [35] has produced the program kaos which allows the 
user to access a variety of algorithms to investigate a differential equa- 
tion from the point of view of modern dynamical systems. 

• There are a number of programs to compute symmetry groups of differ- 
ential equations, including one written by Char [12], and the programs 
SODE and SPDE written by Schwarz [64]. 

• Abelson and Sussman and their group at MIT [1], [2] have combined 
techniques in artificial intelligence to produce software which automat- 
ically analyzes the qualitative features of a differential equation. 

• Wang [72] and Steinberg [68] have developed systems which use sym- 
bolic computation to produce high quality, optimized Fortran code to 
solve differential equations. 
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Figure 1: The trees associated with the vector fields £,. 

• Della- Dora and Tourmer [17] have used the fundamental ideas of Ramis 
[58] to produce a system to analyze linear ordinary differential equa- 
tions. They are now turning their attention to nonlinear systems. 

Point of view. The point of view taken here is to focus on the algorithmic 
aspects of the computation of the vector fields and their brackets and 
to use this information to develop appropriate algorithms which use exact 
symbolic techniques to approximately integrate the trajectories of the dif- 
ferential equation. As will become clear, there are a number of interesting 
points of contact between this approach and the approaches just described. 

In the following sections, we review data structures and algorithms for 
the symbolic computation of the flows of vector fields, and for the symbolic 
approximation of general flows by flows which can be studied symbolically. 

3 Vector fields and the algebra of Cayley trees 

In this section, we describe a data structure which is central to the algo- 
rithms we give for the symbolic computation of senes which approximate 
the solutions of differential equations. The basic idea is to assign trees to 
vector fields as illustrated in Figure 1, and then to impose a multiplication 
on trees which is compatible with the composition of vector fields. 

Consider three vector fields 

E\ = ai D\ + * • • E? =■ b\ D\ + • * • bjqDs * 


£ 3 = c\D\ + . . . c^Ds 

where D, = <9/9x t , and a,, 6 t and c, are smooth functions on R^ Now 
£2 • E\ ~ bj ( Dj ai ) Dj + V] fcjQi Dj D\ 

and £3 * £2 • £1 is equal to 
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Figure 2. The trees associated with Lquation 2, 


4- ^ cj.6 ; a,Dt D ; D, + ^ Ckbj(Dka,)D : D, + ^ c*(D*6 ; )a, D } D,. (2) 

Here the sum is for i, j, k — 1, . . . , N and hence involves 0(N ) differentia- 
tions. It is convenient to keep track of the terms that arise in this way using 
labeled trees: we indicate in Figure 2 the trees that are associated with the 
six sums in this expression. 

An expression such as 

[E Z ,[E 2,£i]] = EiEiE\ - E 3 EiE 2 - E 2 EiE 3 + E x E 2 E 3 (3) 

gives rise in this fashion to 24 trees corresponding to the 24A/^ differenti- 
ations that a naive computation of this expression requires. On the other 
hand. 18 of the trees cancel, saving us from computing 18 N 3 terms. We 
are left with 6 N 3 terms of the form {junk)D^. A careful examination of 
this correspondence between labeled trees and expressions involving the £/ s 
shows that the composition of the vector fields E } ' s, viewed as first order 
differential operators, corresponds to a multiplication on trees. This multi- 
plication is illustrated in Figure 3. It turns out that this construction yields 
an algebra, which we call the algebra of Cayley trees. 

Here is a more precise description for the specialist, following [27] and 
[30]. Let k denote a field of characteristic 0. We say that a finite rooted 
tree is labeled with {Ei,. . . , Em) in case each node, except for the root, is 
assigned a element from this set. Let C'T(E\, ..., Em) denote the set of 
labeled trees and let k{CT{E\, ..., Em)} denote the vector space whose 
basis consists of labeled trees in CT(Ei, ..., Em) ■ Suppose that 1 1 , 
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Figure 3: An example of multiplying two trees. 


t 2 e £T(£j, .... Em) are trees. Let , . . . , s r be the children of the root 
of <i . [f < 2 has n + 1 nodes (counting the root), there are (n + l) r ways to 
attach the r subtrees of <i which have sj , .. . , s r as roots to the tree t 2 by 
making each s; the child of some node of ( 2 . The product t\t 2 is defined 
to be the sum of these (n + l) r trees. It can be shown that this product is 
associative, and that the trivial tree consisting only of the root is a right and 
left unit for this product. In [27], we define a comultiplication on fc{£T(£j, 

. , E,\f )} and show that 

Theorem 3.1 The vector space k {CT {E\ , . . . , Em)} with basis all equiva- 
lence classes of finite rooted trees is a cocommutative graded connected Hop f 
algebra. 

We call this algebra the algebra of Cayley trees generated by the set of 

labeled trees £7(£i Em) The relation between trees and differential 

operators goes back to Cayley [8], [9]. The coalgebra structure on this space 
is very similar to the coalgebra structure defined by Joni and Rota [44]. 
However, the coalgebra structure defined there was defined for individual 
combinatorial objects, rather than for a class of objects such as the family 
of rooted trees. Butcher [6] and [7] has also defined a multiplication on the 
vector space which is dual to the space of trees. This multiplication is closely 
related to the one just defined. 


4 Symbolic evaluation of vector field expressions 

When expressions involving vector fields, such as Lie brackets, are written 
out in coordinates, there is typically a lot of cancellation. Similar cancel- 
lation occurs in expressions involving Poisson brackets and when flows are 
concatenated, as in Campbell-Baker-Hausdorf expansions. In this section, 
we use the algebra of Cayley trees to exploit this cancellation in order to 
compute more efficiently formal expressions involving vector fields. 

This is the set up. Let R denote a ring of functions. In applications, R 
is usually either the ring of polynomial functions, rational functions, or C°° 
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functions. Fix several first order differential operators with coefficients from 

* 

Ej = Y, a 1 D »' a j€«, ( 4 ) 

M=1 

that are defined in terms of a basis of first order differential operators 





ft = 1 , . ■ ■ , N. 


Simplifying any expression in the E : 's using Equation (4) yields a differential 
operator, which we can view as an element of End R. 

We now define a homomorphism from the algebra of expressions in the 
E/s to the algebra of trees. Indeed, the assignment in Figure 1 extends to 
an algebra homomorphism from the free associate algebra in the symbols E : 
to the algebra of Cayley trees k{CT(E\ , . . . , E M )} ■ It is straightforward to 
define a downward-pointing arrow so that the following diagram commutes: 

k<E] , . . . , Em > — * k{CT(E \ , . . , £m)} 

\ 1 ( 5 ) 

End R 


Algorithm 4.1 To rewrite expressions in the first order differential opera- 
tors Ej in terms of the basis 


d 


a 2 

dx jh 0 x 


ft 1 1 i * ■ • — 1 1 , iV, 


compute the composition of the rightward and downward pointing arrows in 
the diagram above. 


In [28], [29] and [33], we show that the algorithm is much more efficient than 
naive substitution, which corresponds to computing the diagonal arrow di- 
rectly. In some common cases, the improvement in efficiency is exponential. 
We have implemented this and related algorithms in Maple, Mathematica 
and Snobol. 

We illustrate this algorithrn by working the example of the last section 
following [31]: consider a higher order derivation of the form 


p = E3E2E1 — E^EiEi — E2E\E$+ E\E2E3. 

Naive simplification requires computing 24 N 3 terms of the form described 
in Table 1. The image of p in the algebra of Cayley trees contains 24 
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No. of terms 

Form of terms 

8 N* ' 

coeff. D M1 

12jV 3 

coeff. 

4jV 3 

coeff. 


Table 1: Naive computation of the differential operator corresponding to p. 



Figure 4: The surviving labeled trees. 


trees, six for each of the four terms of p. For example, the six labeled trees 
corresponding to the first term are given in Figure 2. Eighteen of these 
trees cancel, leaving the six trees in Figure 4. The corresponding differential 
operator is equal to 

4* ^2 a 3 3a 2 2 (^M3^M2 a i tl )D$*\ ~~ T. °3 3 Q 1 2 q 2 )^1» 

and contains 18/V 3 fewer terms^of the form indicated in Table 2 than does 
the naive computation of p. An example of the cancellation of labeled trees 
is given in Figure 5. 

To summarize: we have defined an algebraic structure on families of trees 
which mirrors the algebraic structure of formal expressions in the variables 
Ej , but which alleviates the need for computing intermediate expressions 
which cancel when the noncommuting £/s are expressed in terms of the 
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Figure 5: The term E^E\ £2 contributes the first labeled tree and the term 
E\EiE$ contributes the second, which cancel. 


No. of terms 

Form of terms 

2N 2 3 

coeff. 

12N 3 

coeff. D M D^ t 

AN 3 

coeff. Oji, 


Table 2: Terms in the computation p which cancel. 


commuting D^s. In the following sections, we will see other expressions of 
this simple idea. 

Algorithm 4.1 can be extended in several different directions. 

1. The examples above concern vector fields defined on It is possible 
to work out similar algorithms for vector fields defined on more general 
objects, such as Lie groups. This is important for applications in 
robotics and rigid body dynamics. For example, the group G — 50(3) 
is the appropriate configuration space for a rotating rigid body. To be 
more specific, assume the vector fields are of the form 

E, = £ 

M=1 

where the Y 4 are left-invariant vector fields on G and the are smooth 
functions on the group. In this case, the Cayley algebras are generated 
by ordered , labeled trees {16]. Roughly speaking, the trees are ordered 
since the vector fields no longer commute. 

2. The natural action of differential operators on functions turns the ring 

of functions R into a module. In this same way, the trees have a natural 
action on R , as indicated in the Diagram 5. It turns out [34] that this 
gives R the structure of /f/t-bialgebra, as introduced by Nichols [55] 
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and [56]. These types of algebras are closely related to differential 
algebras. 

3. It is a basic fact that the local properties of the nonlinear system ( 1) are 
determined by the algebraic properties of the higher order iterated Lie 
brackets; see, for example, [40]. Unfortunately, due to intermediate 
expansion swell, it is often difficult to compute these using current 
computer algebra systems. It turns out that higher order Lie brackets 
not only involve the cancellation of all terms above the first order but 
also the cancellation of some of the first order terms. Algorithm 4.1 
can be used so that the terms arising in these first order cancellations 
need not be computed. 

5 Exponentials, Lie brackets, and nil flows 

Some differential equations have the property that their flows can be inte- 
grated symbolically in closed form. For example, this holds for differential 
equations of the form (1), if E\ and E? are homogeneous in the appropriate 
sense and generate a graded, niipotent Lie algebra. In this section, we give 
an algorithm which, given an appropriate niipotent Lie algebra, returns vec- 
tor fields on with polynomial coefficients which generate the Lie algebra. 
This leads to an interesting class of differential equations whose flows can 
be integrated symbolically in closed form. At the end of the section, we look 
at several applications of this algorithm. 

By the third fundamental theorem of Lie [63], we know that a niipotent 
Lie algebra arises from some Lie algebra of vector fields. What is not obvious 
is how to construct such vector fields. Niipotent Lie algebras of vector fields 
have been used as an important tool in partial differential equations by 
Folland and Stein [21], Rothschild and Stein [60], and Rockland [61]; and in 
control theory by Krener [48], Hermes [41], [42], and Crouch [15], [14]. We 
will see below how they are also a useful tool in developing symbolic-numeric 
algorithms to integrate flows. 

Our goal is to describe a natural representation of niipotent Lie algebras 
on vector fields on Euclidean space with polynomial coefficients. To define 
this represenation, we define a basts of Hall trees on generators E\ and £2 
recursively as follows: 

1. basis elements consist of rooted, binary trees, with all nodes, except 
the root, labeled with E\ y £ 2 , £ 3 , satisfying 
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(a) all right children are leaves and labeled with either E\ or £ 2 

(b) the sequence formed by the labels of the right leaves at increasing 
distance from the root is nonincreasing 

2. the two rooted trees consisting of a root and a single (left) child labeled 

for i = 1,2 are in the basis and of length 1 

3. if we have defined basis elements t x of length 1 , . . ., r - 1 , they are 
simply ordered so that t t < t } in case the length of t t is less than the 
length of t y 

4. a tree consisting of a root with a single (left) child labeled E t is in the 
basis provided that E t ' s left child is in the basis and of lower order. 

To each such tree corresponds an element E x in the Hall basis [39] of the 
free, nilpotent Lie algebra on two generators. Figure 6 illustrates how the 
basis of Hall trees leads naturally to a representation of the free nilpotent 
Lie algebra generated by E x and £ 2 on the space of vector fields on R N 
with polynomial coefficients. See [22] for further details of the following 
algorithm: 

Algorithm 5.1 Fix r > 1 and say that the free, nilpotent Lie algebra of 
two generators of rank r has dimension jV. Let E t , . . . , Ejj denote the Hall 
basis. Then the vector fields on R.^ defined as in Figure 6 satisfy 

1. the Lie algebra they generate is isomorphic to the free, nilpotent Lie 
algebra on two generators of rank r 

2. any trajectory of the nonlinear system 

x(t) = E,(*(0) + u(0 E 2 (i( 0), *(0) = x° € 

can be computed in closed form in terms of quadratures of the function 
t - U (0 

3- there exist constants o\, . . qjj such that 



exp(E 2 )e\p(Ei) = exp 
and the a, can be computed by solving a lower triangular linear system. 


13 




Figure 6: The vector fields arising from the basis of Hall trees. The first 
tree is sent to D \ , while the sum of the trees is sent to £>2— *1 £>3 D4 

+ xix 2 D s - - \x\x 2D 7 -\x\x\Dg. 


i 
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As a note, tliis algorithm was found only after several months of experimen- 
tation with Maple and required the careful study of the free, nilpotent Lie 
algebra on two generators of rank 5, which is 23 dimensional. 

We conclude this section with some remarks describing several applica- 
tions of this algorithm and related work. 

1. Locally approximating a nonlinear system by an explicitly integrable 
nilpotent one yields a number of integration algorithms, which we re- 
fer to as pieceiLise nilpotent integration algorithms (PWNI). The basic 
idea [23], [21] is to approximate locally a nonlinear system at a given 
point by an explicitly integrable nilpotent one in which the computa- 
tions can be done symbolically in closed form, and to patch together 
the various nilpotent approximations at nearby points using a stan- 
dard numerical algorithm. Preliminary work indicates that this leads 
to symbolic-numeric algorithms for the path planning problem and ef- 
ficient algorithms to integrate neighborhoods of trajectories around a 
given fixed reference trajectory. 

2. Algorithm 5.1 also provides an efficient means for computing the con- 
catenation of flows. Write z(t) = exp Ex 0 for the flow of the nonlinear 
system 

*(<)=£(*(*)), *( 0 ) = *°, ( 6 ) 

The Campbell-Baker-FIausdorff formula expresses the concatenation 
of two flows as a single flow; 

exp(tE 2 )exp(tE\) = exp (tE 2 +tE\ + l/2f 2 [E 2l E\] 


+ l/12p2> Ei], Ei) - 1/12[[E 2 , Ei] t E 2 ) + ■ ■ •)■ 


Consider the equation 


(7) 


N 

exp (E 2 )exp(E l ) = ^c t E lt 

i=l 

where E t are the vector Ifields produced by the algorithm. Since all 
flows of the vector fields E\ and E 2 are explicitly integrable in closed 
form, this reduces the computation of the c, to the solution of a linear 
system, which is lower triangular. 

3. We can also use Algorithm 5.1 to derive a class of numerical integra- 
tors, which are sometimes known as splitting methods. For example, 
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suppose that £, and £ 2 are separately integrable in closed form, but 
E\ + £2 is not. Then using the algorithm, we can compute constants 
r, such that 

exp(r 7 £i) exp(r 6 £ 2 ) • exp(r 5 £i) • exp(r 4 £ 2 ) 

■ exp(r 3 £i) • exp(r 2 £ 2 ) • exp(n £i)- = exp(£i + £ 2 ) + 0(t 5 ) (8) 

This formula yields a numerical integrator for our original nonlinear 
system (1) (with u = 1) This algorithm is used in accelerator physics 

[65]. 

4. Recently, Strichartz [69] has shown that the solution of the initial value 
problem 

x(t) = E(t,x(t)), x(0) = x° 

can be written as 

z(t) = exp(G(0) x °> 

with 

OO 

r = l Q 

‘ J [‘ ' ' (£*( s <y(l) 1 ^( s <r(2)l * ’ *]^( S ^(r)] 

and where <j ranges over the symmetric group on r symbols, the inte- 
gration region is a simplex in R r , £(s) denotes the vector field E(s , ■), 
and exp is defined as in Equation 6. Formulas of this type date back 
to Chen [13] and appear to be related to Algorithm 5.1. 

5. F. Bergeron, N. Bergeron, and A. M. Garsia have also exploited a rela- 
tion between trees and polynomials in their study of free Lie algebras; 
see [4] and the references cited there. 

6 Taylor series and intrinsic integrators 

Although nonlinear systems often conserve quantities such as energy or an- 
gular momentum, most numerical integrators do not. Similarly, nonlinear 
systems typically evolve on some underlying geometric space, such as a Lie 
group or homogeneous space, but most numerical integrators do not remain 
in such a space. 

Recently, there has been a flurry of activity related to numerical in- 
tegrators preserving the symplectic structure, sparked off by the work of 
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Channell [10] and Scovel [11], and based upon earlier work of deVogeieaere 
[19]. These types of numerical schemes have found important applications in 
the long term study of orbits in accelerator physics and in other areas [65]. 
The derivation of these integrators typically involves the symbolic computa- 
tion of Taylor series and generating series for the symplectic transformation 
which is the update for the numerical integrator. 

Consider a numerical integrator for a differential equation 

i(<) = E(*(0), *(0 ) = p€Af 

evolving on a space Xf . Call a numerical integrator intrinsic in case z n £ Xf 
implies r n +i € Xf , for n > 1, where r n is the approximation to the trajectory 
at time t n . One means of deriving intrinsic numerical integrators is 
to rnimic the derivation of standard numerical integrators, but to impose 
additional constraints on the scheme to satisfy the added condition that the 
points x n remain in the space A/. This typically involves the careful study 
of the Taylor series of the solution. 

This can be done by using the Cayley algebra of trees, as briefly indicated 
in [32]. As an illustration of this, we consider intrinsic Runge-Kutta type 
algorithms evolving on a Lie group G, following [16]. Let g denote its Lie 
algebra, and let Y\ , . Ytf denote a basis of g . We give an algorithm to 
approximate solutions to differential equations evolving on G of the form: 

i(t) = E(z{t)), x(Q ) = p€G, 

where 

X 

* 4=1 

and the = a M are analytic functions on G. Let exp(f£*) • z denote the 
solution z(t) at time t. The algorithms depends upon constants c, and c l; , 
for 1 = 1- • • ■ . £ and j < i. For fixed constants, define the following elements 
of the Lie algebra g 

X 

El = <? 

*4 = 1 

X 

E 2 = ^2 a u (exp{hc 2 iEi) ■ € g 

* 4=1 

X 

^3 = yi QM (exp(/ic 3 2 E 2 ) exp(ftc^i£i) fa... 

M=1 
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Figure 7 Trees associated with third order terms in a Taylor series. 

These arise by “freezing the coefficients” of E at various points along the 
flow of E . 

Algorithm 6.1 Given an initial point xq on the group, define 
x n + 1 = exp hckEk • exp hc } E\X n , 

for n > 0. 

Notice that if we assume the exponential exp(/i£,) maps the Lie algebra 
to the Lie group exactly, then this algorithm is intrinsic. For a group such as 
G — 50(3), there are classical closed form expressions for the exponential 
map and Algorithm 6.1 yields an intrinsic integrator. Notice also that if G 
is the abelian group R^, then the algorithm becomes the classical Runge- 
Kutta algorithm. 

The first step is to derive the equations that the coefficients c t and c tJ 
must satisfy in order for the algorithm to yield an rth order numerical in- 
tegrator. This can easily be done using the Cayley algebra of ordered trees 
[16]. The trees are ordered since the vector fields do not commute. 

Assume for the moment that G = R^ and consider the terms in the 
Taylor series 

z(t+h)-z(t) = hz(t) + ^i(t) + ^x {u,) {t)+ ■ 

h 2 

= hE+—DEE 

h 3 f 2 \ 

+— (i DE DE E + D 2 E{E , £)) + • • 

Notice that there is a natural correspondence between trees and terms in 
the series. For example, the /i 3 /3! terms are associated with the two labeled 
trees in Figure 7. This observation goes back to at least Cayley [8], [9]. 

We now generalize this to a Lie group following [16]. Recall that Dia- 
gram 5 induces an action of trees on the ring of analytic functions on G. 
Using this action, we can now state the 
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Lenuna G.l Lei a denote (he tree consisting of a root with a single child 
labeled F Then for any analytic function f on the group and for sufficiently 
small t, 

f(exp(tF) ■ x) = exp(Ia) • / U ■ 

Notice that if G is Euclidean space, and if the functions / are the coordinate 
functions xj x„, then this becomes the familar Taylor series. 

Using this lemma, it is now easy to compute the equations that the 
coefficients c, and c tJ must satisfy. In spirit, this is similar to Butcher’s use 
of trees to analyze higher order Runge-Kutta algorithms in Euclidean space 

[ 6 ], [ 71 - 
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